#include "CompBioTools.hh"
#include "MerHashMap.hh"


int OptMerSubset =0; // -1 means do them all. 
//string OptMerSubset ="ALL"; // 0 means do them all. 
vector<string> OptFiles;
int OptMaxMer = 3;
string OptOutFileName = "";
string OptFON = "";
string OptDistOutName;
 
//---------------------------------------------------
// * PrintHelp
//---------------------------------------------------
//
void PrintHelp(){
printf("\n
    atlas-abundant-mer -F FON -R {MaxMer} -s {Subset} -o {mer output file} -d {dist out} 

    If no FON is given reads from standard input.  

    The output is the set of mers that occur more than MaxMer times. 
    Right now the mer size is fixed at 32-mer.  In the future this will
    probably be optional. 


    Subset has been changed to be a number 0-63.  There are 64 batches. 
    It is hoped that any given batch will be < 100 million uncollapsed 
    (50 million collapsed) mers in order to fit in the 1GB limit at 
    10 bytes/uncollapsed mer.  If the sequence contains the worst case 
    of 2.75*10^9*2 = 5.5*10^9 distinct uncollapsed mers, divided by 
    64 = 86 million mers/batch.  It's going to be close. 

    Entering a subset of 999 results in all of the sequence being hashed. 


    Written by:  James Durbin (kdurbin@bcm.tmc.edu)
        
    Version:  1.00

");
}

//---------------------------------------------------
// * SetupOptions
//---------------------------------------------------
//
void SetupOptions(int argc, char**argv)
{
  // Print it for the record
  for(int i = 0;i< argc;i++) cerr<<argv[i]<<" ";
  cerr<<endl;

  // Create parser object and define valid options
  // Options with a : require a parameter. 
  GetOption Parser(argc,argv,"R:d:s:o:F:h");

  // Handle the options...
  char theOption;
  while((theOption = Parser.NextOption()) != OPTION_EOF){
    switch(theOption){
	 case 'F': OptFON = Parser.GetStringArg();break;
	 case 'R': OptMaxMer = Parser.GetIntArg();break;
	 case 's': OptMerSubset = Parser.GetIntArg();break;
	 case 'o': OptOutFileName = Parser.GetStringArg();break;
	 case 'd': OptDistOutName = Parser.GetStringArg();break;
    default: PrintHelp(); exit(1); break;
    }
  }

  // Add query files on command line to list of query files.
  string temp = Parser.GetFreeArg();
  while(temp != ""){
	 OptFiles.push_back(temp);
	 temp = Parser.GetFreeArg();
  }
}

int main(int argc,char** argv)
{

  SetupOptions(argc,argv);

  Mer2IntMap Mers;

  if (OptFON !=""){
	 ReadFOF(OptFON,OptFiles);
  }

  // If no files were specified, then read from the input stream
  if (OptFiles.size() == 0){
	 if (OptMerSubset == -1){
		Mers.Add(cin);
	 }else{
		Mers.Add(cin,OptMerSubset);
	 }
  }else{
	 
	 for(int i = 0;i< OptFiles.size();i++){		
		izipstream fin;
		fin.open(OptFiles[i].c_str());
		if (fin.fail()){
		  cerr<<"Couldn't open "<<OptFiles[i]<<endl;
		}
		cerr<<"Adding "<<OptFiles[i]<<" to hash map..."<<flush;
		
		if (OptMerSubset == 999){
		  Mers.Add(fin);
		}else{
		  Mers.Add(fin,OptMerSubset);
		}
		cerr<<"\t"<<Mers.size()<<" distinct mers.\t"<<Mers.mTotalMersInHash<<" total mers."<<endl;
		fin.close();
	 }
  }

  ofstream fout;
  OpenOrDie(OptOutFileName,fout);
  
  // Compute stats 
  
  cout<<"Computing statistics..."<<endl;
  int MaxMer = 0;
  long TotalMerCount = 0;
  Mer2IntMap::iterator it;

  vector<int> MerDist;
  int MaxMerCount = 200000;
  int NumOutOfRange = 0;
  MerDist.resize(MaxMerCount);
  for(int i = 0;i< MaxMerCount;i++){
	 MerDist[i] = 0;
  }

  for(it = Mers.begin();it!=Mers.end();it++){

	 int merCount = (*it).second;

	 if (merCount < MaxMerCount){
		MerDist[merCount]++;
	 }else{
		NumOutOfRange++;
	 }

	 if (merCount > MaxMer) MaxMer = merCount;
	 
	 TotalMerCount+=merCount;
	 
	 if (merCount > OptMaxMer){
		Sequence Seq;
		Mer TempMer = (*it).first;
		TempMer.ToSequence(Seq);
		fout<<Seq<<"\t"<<merCount<<endl;
	 }		  
  }
  fout.close();


  double MeanMerFreq = (double)TotalMerCount/(double)Mers.size();

  double DeviationSum = 0;
  for(it = Mers.begin();it!=Mers.end();it++){

	 int merCount = (*it).second;
	 double deviation = MeanMerFreq - merCount;
	 double devsq = deviation*deviation;
	 DeviationSum +=devsq;
  }

  double Variance = DeviationSum/(double)Mers.size();
  double stdev = sqrt(Variance);

  ofstream dout;
  dout.open(OptDistOutName.c_str());
  
  dout<<OptMerSubset<<"\t"<<Mers.size()<<"\t"<<Mers.mTotalMersInHash;
  dout<<"\t"<<MeanMerFreq<<"\t"<<stdev<<"\t"<<MaxMer<<endl;

  for(int i = 0;i< MaxMerCount;i++){
	 dout<<i<<"\t"<<MerDist[i]<<endl;
  }
  dout<<">="<<MaxMerCount<<"\t"<<NumOutOfRange<<endl;

}


